#include "jumppaths.hpp"

jumppaths::~jumppaths(){
  delete[] tmap;
}

jumppaths::jumppaths(const std::vector<matrix<int> > &ops, const vector<std::complex<double> > &psi){
  Psi=psi;
  Ops=ops;
  int n=Ops.size();
  tmap=new transition_map[n];
};
bool jumppaths::is_orthogonal(vector<std::complex<double> > vec1,vector<std::complex<double> >vec2){
  bool ret;
  //If innerproduct is smaller than 1e-10 then vectors are considered orthogonal
  double thres=1e-10;
  if(abs(inner_prod(vec1,vec2))<thres){
    ret=true;
  }
  else{
    ret=false;
  }
  return ret;
}

bool jumppaths::is_zero_vector(vector<std::complex<double> > vec){
  bool ret;
  //If sum of real and imag parts of all elements is < thres, then zero vector
  double thres=1e-10;
  double sum=0;
  vector<std::complex<double> >::iterator it;
  for(it=vec.begin();it!=vec.end();it++){
    sum+=std::real(*it)+std::imag(*it);
  }
  if(sum<thres){
    ret=true;
  }
  else{
    ret=false;
  }
  return ret;
}



void jumppaths::calculate_transition_map(){
  //Initial state is indexed with zero;
  int state_num=0;
  int num_of_channels=Ops.size();
  std::vector<vector<std::complex<double> > > states_old,states_new;
  //Initialization
  states_old.assign(Ops.size(),Psi);
  states_new.assign(Ops.size(),Psi);
  std::vector<matrix<int> >::iterator op_it;
  std::vector<vector<std::complex<double> > >::iterator vec_it_new=states_new.begin();
  std::vector<vector<std::complex<double> > >::iterator vec_it_old=states_old.begin();
  
  if(states_new.size()!=Ops.size()||states_old.size()!=Ops.size()){
    std::cerr<<"Error"<<std::endl;
    exit(1);
  }
  
  for(op_it=Ops.begin();op_it!=Ops.end();op_it++,vec_it_new++,vec_it_old++){
    *vec_it_new = (*op_it) * (*vec_it_old);
  }


};
